Measurement exploratory data analysis¶

In [ ]:
DEVICES = [
    'beaglebone-fan',
    'beaglebone-compressor',
    'beaglebone-pump',
    'beaglebone-refrigerator',
    'mafaulda-a',
    'mafaulda-b'
]

DEVICE = DEVICES[3]
T_WAVEFORM = 5  # (1 = MaufaulDa, 5 = others)
T_SEC = T_WAVEFORM
NFFT = 2**10   # (2**10 = MaufaulDa, 2**14 = others)
F_LIMIT = None # (3000 = MaufaulDa, None = others)
In [ ]:
import os
import pandas as pd
import numpy as np
import matplotlib.pylab as plt

from tabulate import tabulate
from IPython.display import Markdown, HTML
from tqdm.notebook import tqdm

from typing import List, Tuple
from scipy.signal import find_peaks, butter, lfilter
from tsfel.feature_extraction.features import fundamental_frequency

from zipfile import ZipFile
import sys
sys.path.append('../')
from vibrodiagnostics import mafaulda, selection, discovery, models


def beaglebone_measurement(filename: str, fs: int) -> Tuple[str, pd.DataFrame]:
    g = 9.81
    milivolts = 1800
    resolution = 2**12
    columns = ['x', 'y', 'z']
    ts = pd.read_csv(filename, delimiter='\t', index_col=False, header=None, names=columns)
        
    # Calculate amplitude in m/s^2 Beaglebone Black ADC and ADXL335 resolution (VIN 1.8V, 12bits)
    for dim in columns:
        ts[dim] = ts[dim] * (milivolts / resolution)  # ADC to mV
        ts[dim] = (ts[dim] / 180) * g                 # mV to m/s^2 (180 mV/g)
        ts[dim] -= ts[dim].mean()

    ts['t'] = ts.index * (1 / fs)
    ts.set_index('t', inplace=True)
    return (os.path.basename(filename), ts, fs, ts.columns)  # last is feature columns


def beaglebone_dataset(filenames: List[str], fs: int) -> List[Tuple[str, pd.DataFrame]]:
    dataset = []
    for filename in filenames:
        name, ts, fs, cols = beaglebone_measurement(filename, fs)
        dataset.append((name, ts))
    return dataset


def lowpass_filter(data, cutpoint, fs, order=5):
    b, a = butter(order, cutpoint, fs=fs, btype='lowpass')
    y = lfilter(b, a, data)
    return y


def mafaulda_dataset(
        place=['ax', 'ay', 'az'],
        features_path =  '../../datasets/features_data/',
        mafaulda_path='../../datasets/MAFAULDA.zip',
        rpm=2500,
        lowpass_hz=10000):

    metadata_filename = os.path.join(features_path, selection.MAFAULDA_METADATA)
    faults = {
        'normal': 'normal',
        'imbalance': 'imbalance',
        'horizontal-misalignment': 'misalignment',
        'vertical-misalignment': 'misalignment',
        'overhang-cage_fault': 'cage fault',
        'underhang-cage_fault': 'cage fault',
        'underhang-ball_fault': 'ball fault',
        'overhang-ball_fault': 'ball fault',
        'overhang-outer_race': 'outer race fault',
        'underhang-outer_race': 'outer race fault',
    }

    metadata = pd.read_csv(metadata_filename, index_col='filename')
    metadata.reset_index(inplace=True)
    metadata = models.fault_labeling(metadata, faults)
    files = pd.DataFrame()
    # Worst severity and mid rpm
    for name, group in metadata[metadata['rpm'] >= rpm].groupby(by='fault', observed=False):
        files = pd.concat([
            files,
            group[
                group['severity_level'] == group['severity_level'].max()
            ].sort_values(by='rpm', ascending=True).head(1)
        ])
    ordering = {
        'normal': 0,
        'misalignment': 1,
        'imbalance': 2,
        'cage fault': 3,
        'ball fault': 4,
        'outer race fault': 5,
    }
    source = ZipFile(mafaulda_path)
    dataset = len(files) * [0]
    for index, file in files.iterrows():
        ts = mafaulda.csv_import(source, file['filename'])
        ts = ts[place]
        ts.columns = ts.columns.str.extract(r'(\w)$')[0]
        for axis in ts.columns:
            ts[axis] = lowpass_filter(ts[axis], lowpass_hz, file['fs'])
        pos = ordering[file['fault']]
        dataset[pos] = ((file['fault'] + ' (' + file['filename'] +')', ts))

    return dataset

Load dataset

In [ ]:
if DEVICE == 'beaglebone-fan':
    Fs = 2500
    path = '../../inspections/fan/'
    files = [
        '1_still.tsv', '2_still.tsv', '3_still.tsv',
        '1_up.tsv', '2_up.tsv', '3_up.tsv',
        '1_down.tsv', '2_down.tsv', '3_down.tsv'
    ]
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-compressor':
    Fs = 2500
    path = '../../inspections/datacentres/shc3/'
    files = [
        'k3_1.tsv', 'k3_2.tsv', 'k3_3.tsv', 'k3_4.tsv',
        'k5_1.tsv', 'k5_2.tsv', 'k5_3.tsv', 'k5_4.tsv'
    ]
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-pump':
    Fs = 2500
    path = '../../inspections/bvs/'
    files = [
        'bvs_1_hore.tsv', 'bvs_2_hore.tsv' 
        #, 'bvs_3_motor.tsv', 'bvs_4_motor.tsv'
    ]    
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-refrigerator':
    Fs = 2500
    path = '../../inspections/home-refrigerator/'
    files = [
        'ch1.tsv', 'ch2.tsv', 'ch3.tsv', 'ch4.tsv', 'ch5.tsv'
    ]    
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'mafaulda-a':
    Fs = mafaulda.FS_HZ
    DATASET = mafaulda_dataset(place=['ax', 'ay', 'az'])

elif DEVICE == 'mafaulda-b':
    Fs = mafaulda.FS_HZ
    DATASET = mafaulda_dataset(place=['bx', 'by', 'bz'])
In [ ]:
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts.info()
    print()

ch1.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 25400 entries, 0.0 to 10.159600000000001
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       25400 non-null  float64
 1   y       25400 non-null  float64
 2   z       25400 non-null  float64
dtypes: float64(3)
memory usage: 793.8 KB

ch2.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 25400 entries, 0.0 to 10.159600000000001
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       25400 non-null  float64
 1   y       25400 non-null  float64
 2   z       25400 non-null  float64
dtypes: float64(3)
memory usage: 793.8 KB

ch3.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 25400 entries, 0.0 to 10.159600000000001
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       25400 non-null  float64
 1   y       25400 non-null  float64
 2   z       25400 non-null  float64
dtypes: float64(3)
memory usage: 793.8 KB

ch4.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 25400 entries, 0.0 to 10.159600000000001
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       25400 non-null  float64
 1   y       25400 non-null  float64
 2   z       25400 non-null  float64
dtypes: float64(3)
memory usage: 793.8 KB

ch5.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 25400 entries, 0.0 to 10.159600000000001
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       25400 non-null  float64
 1   y       25400 non-null  float64
 2   z       25400 non-null  float64
dtypes: float64(3)
memory usage: 793.8 KB

In [ ]:
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    display(tabulate(ts.describe(), headers='keys', tablefmt='html'))
    ts.boxplot(grid=True)
    plt.show()

ch1.tsv

x y z
count25400 25400 25400
mean -3.87917e-15 2.07876e-15 -5.41299e-16
std 0.0684914 0.0667617 0.0383097
min -0.2103 -0.219397 -0.145902
25% -0.0665985 -0.0517456 -0.026151
50% 0.00525207 -0.00384523 -0.00220078
75% 0.0531525 0.0440552 0.0217494
max 0.196854 0.211707 0.189401
No description has been provided for this image

ch2.tsv

x y z
count25400 25400 25400
mean -7.03521e-15 3.25031e-15 4.41879e-15
std 0.0658363 0.0677933 0.0384159
min -0.169978 -0.19035 -0.153563
25% -0.0502275 -0.0466491 -0.0338122
50% -0.00232713 0.00125126 -0.00986201
75% 0.0455733 0.0491516 0.0140882
max 0.237175 0.192853 0.157789
No description has been provided for this image

ch3.tsv

x y z
count25400 25400 25400
mean -5.84659e-16 -2.04491e-15 1.47004e-15
std 0.0751321 0.0744217 0.0376824
min -0.239868 -0.238068 -0.14983
25% -0.0722164 -0.0704164 -0.0300792
50% -0.000365853 0.00143418 -0.00612899
75% 0.0714847 0.0493346 0.0178212
max 0.215186 0.216986 0.161522
No description has been provided for this image

ch4.tsv

x y z
count25400 25400 25400
mean -2.22031e-15 1.16932e-16 -7.02318e-15
std 0.0744901 0.0851 0.0369085
min -0.217635 -0.280152 -0.139479
25% -0.0499833 -0.0646005 -0.0197278
50% 0.0218673 0.00725012 0.0042224
75% 0.0458175 0.0791007 0.0281726
max 0.213469 0.222802 0.147924
No description has been provided for this image

ch5.tsv

x y z
count25400 25400 25400
mean 6.38929e-16 9.97557e-15 -9.22587e-16
std 0.0675674 0.0792405 0.0374805
min -0.216875 -0.254568 -0.132351
25% -0.0492233 -0.0629664 -0.0365504
50% -0.00132292 -0.0030909 0.0113499
75% 0.0465775 0.0567846 0.0353001
max 0.190279 0.248386 0.155051
No description has been provided for this image

Statistical tests

  • Normality test: Kolmogorov–Smirnov test
  • Normality visual test: Quantile-quantile plot on chosen recording
  • Stationarity test: Augmented Dickey–Fuller test
  • Stationarity visual test: Autocorrelation plot
In [ ]:
from statsmodels.tsa.stattools import adfuller
from statsmodels.api import qqplot
from scipy.stats import kstest

normality_tests = []
for name, ts in DATASET:
    for x in ts.columns:
        p_value = kstest(ts[x], 'norm').pvalue
        test = {'name': name, 'axis': x, 'p-value': p_value, 'not-normal': p_value < 0.05}
        normality_tests.append(test)

normality_tests = pd.DataFrame.from_records(normality_tests)
print(normality_tests.value_counts('not-normal'))
normality_tests.describe()
not-normal
True    15
Name: count, dtype: int64
Out[ ]:
p-value
count 15.0
mean 0.0
std 0.0
min 0.0
25% 0.0
50% 0.0
75% 0.0
max 0.0
In [ ]:
name, ts = DATASET[0]
fig, ax = plt.subplots(1, len(ts.columns), figsize=(10, 4))
for i, x in enumerate(ts.columns):
    qqplot(ts[x], line='45', ax=ax[i], marker='.', alpha=0.5)
    ax[i].set_title(f'Axis: {x}')

plt.tight_layout()
print(name)
plt.show()
ch1.tsv
No description has been provided for this image
In [ ]:
stationarity_tests = []
for name, ts in tqdm(DATASET):
    for x in ts.columns:
        result = adfuller(ts[x].loc[T_WAVEFORM:T_WAVEFORM+1])
        p_value = result[1]
        test = {
            'name': name,
            'axis': x,
            'statistic': result[0],
            'p-value': p_value,
            'stationary': p_value < 0.001
        }
        stationarity_tests.append(test)

stationarity_tests = pd.DataFrame.from_records(stationarity_tests)
print(stationarity_tests.value_counts('stationary'))
stationarity_tests['p-value'].describe()
  0%|          | 0/5 [00:00<?, ?it/s]
stationary
True    15
Name: count, dtype: int64
Out[ ]:
count    1.500000e+01
mean     3.856455e-17
std      1.493597e-16
min      0.000000e+00
25%      0.000000e+00
50%      2.330396e-30
75%      1.487401e-23
max      5.784678e-16
Name: p-value, dtype: float64
In [ ]:
name, ts = DATASET[0]
fig, ax = plt.subplots(1, len(ts.columns), figsize=(10, 4))
for i, x in enumerate(ts.columns):
    ax[i].acorr(ts[x], maxlags=50)
    ax[i].set_title(f'Axis: {x}')

plt.tight_layout()
print(name)
plt.show()
ch1.tsv
No description has been provided for this image

Time domain histogram

In [ ]:
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    axis = ts.columns
    ax = ts[axis].hist(figsize=(15, 3), grid=True, bins=100, layout=(1, 3), edgecolor='black', linewidth=0.5)
    plt.show()

ch1.tsv

No description has been provided for this image

ch2.tsv

No description has been provided for this image

ch3.tsv

No description has been provided for this image

ch4.tsv

No description has been provided for this image

ch5.tsv

No description has been provided for this image

Time domain waveform

In [ ]:
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    axis = ts.columns
    
    ax = ts[axis].plot(figsize=(20, 8), grid=True, subplots=True)
    for i, axname in enumerate(axis):
        ax[i].set_xlabel('Time [s]')
        ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
    plt.show()               # plt.savefig('waveform.png')

ch1.tsv

No description has been provided for this image

ch2.tsv

No description has been provided for this image

ch3.tsv

No description has been provided for this image

ch4.tsv

No description has been provided for this image

ch5.tsv

No description has been provided for this image

Time domain waveform zoom detail

In [ ]:
for name, ts in DATASET:
    axis = ts.columns
    display(Markdown(f'**{name}**'))
    ax = (ts[axis].iloc[int(T_WAVEFORM*Fs):int(T_WAVEFORM*Fs)+Fs]
                  .plot(figsize=(20, 10), grid=True, subplots=True))
    
    for i, axname in enumerate(axis):
        ax[i].set_xlabel('Time [s]')
        ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
        plt.show()      # plt.savefig('waveform_zoom.png')

ch1.tsv

No description has been provided for this image

ch2.tsv

No description has been provided for this image

ch3.tsv

No description has been provided for this image

ch4.tsv

No description has been provided for this image

ch5.tsv

No description has been provided for this image

Time domain waveform zoom - faults side by side

In [ ]:
fig, ax = plt.subplots(len(DATASET), 3, figsize=(12, 15), sharex=True)

for idx, df in enumerate(DATASET):
    name, ts = df
    columns = ts.columns
    ax[idx][1].set_title(name)
    ax[idx][0].set_ylabel('Amplitude [m/s^2]')

    for pos, axis in enumerate(columns):
        data = ts[axis].loc[T_WAVEFORM:T_WAVEFORM+0.3]
        ax[idx][pos].plot(data.index, data, linewidth=1, color='darkblue')
        ax[idx][pos].grid()
    

plt.tight_layout()
plt.show()
No description has been provided for this image
In [ ]:
def spectogram(x, debug=True):
    fig, ax = plt.subplots(figsize=(15, 4))
    cmap = plt.get_cmap('inferno')
    pxx, freqs, t, im = plt.specgram(
        x, NFFT=NFFT, Fs=Fs,
        detrend='mean',
        mode='magnitude', scale='dB',
        cmap=cmap, vmin=-60
    )
    fig.colorbar(im, aspect=20, pad=0.04)
    ax.set_xlabel('Time [s]')
    ax.set_ylabel('Frequency [Hz]')
    mafaulda.resolution_calc(Fs, NFFT)
    return freqs, pxx


def window_idx(t):
    return (Fs * t) // NFFT + 1


def spectrum_slice(freqs, Pxx, t):
    fig, ax = plt.subplots(2, 1, figsize=(20, 8))
    n = window_idx(t)

    dB = 20 * np.log10(Pxx.T[n] / 0.000001)
    ax[0].plot(freqs, dB)      # 1 dB = 1 um/s^2
    ax[0].grid(True)
    ax[0].set_xlabel('Frequency [Hz]')
    ax[0].set_ylabel('Amplitude [dB]')
    
    ax[1].plot(freqs, Pxx.T[n])
    ax[1].grid(True)
    ax[1].set_xlabel('Frequency [Hz]')
    ax[1].set_ylabel('Amplitude [m/s^2]')
    return n


def get_max_frequency(freqs, Pxx, i):
    max_freq = freqs[np.argmax(Pxx.T[i])]
    return max_freq


def get_peaks(freqs, Pxx, i, top=5):
    amplitudes = Pxx.T[i]
    peaks, _ = find_peaks(amplitudes, distance=3)

    fundamental = get_max_frequency(freqs, Pxx, i)
    f_top = freqs[peaks[np.argsort(amplitudes[peaks])]][::-top]
    y_top = np.sort(amplitudes[peaks])[::-top]

    return pd.DataFrame({
        'f': f_top,
        'y': y_top,
        '1x': f_top / fundamental 
    })


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter(order, [lowcut, highcut], fs=fs, btype='band')
    y = lfilter(b, a, data)
    return y


def get_spectrograms(DATASET: List[pd.DataFrame], axis: str) -> list:
    spectrograms = []

    for name, ts in DATASET:
        base_freq = fundamental_frequency(ts[axis], Fs)
        display(Markdown(f'**{name}** *({axis.upper()} axis, Fundamental = {base_freq:.4f} Hz)*'))
        
        freqs, Pxx = spectogram(ts[axis])
        spectrograms.append((name, freqs, Pxx))
        plt.show()          # plt.savefig(f'x_axis_fft_{NFFT}.png')
    
    return spectrograms


def show_spectrogram_detail(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
        i_window = spectrum_slice(freqs, Pxx, t)
        plt.show()           #plt.savefig(f'x_axis_fft_{NFFT}_at_{T_SEC}s.png')


def show_mms_peaks(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
    
        i_window = window_idx(t)
        peaks = discovery.mms_peak_finder(Pxx.T[i_window])
        
        fig, ax = plt.subplots(1, 1, figsize=(15, 3))
        ax.grid(True)
        ax.plot(freqs, Pxx.T[i_window])
        ax.scatter(freqs[peaks], Pxx.T[i_window][peaks], marker='^', color='red')
        ax.set_xlabel('Frequency [Hz]')
        
        plt.show()


def show_harmonic_series(spectrograms: list, axis: str, t: float):
    # https://stackoverflow.com/questions/1982770/changing-the-color-of-an-axis
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
    
        i_window = window_idx(t)
        h_series = discovery.harmonic_series_detection(freqs, Pxx.T[i_window], Fs, NFFT)
    
        # Find best (sum of harmonics' amplitudes in the largest)
        max_harmonic_amp_idx = np.argmax([
            sum([h[1] for h in s]) / len(s)
            for s in h_series
        ])
        best_harmonic_series = pd.DataFrame(
            h_series[max_harmonic_amp_idx],
            columns=['Frequency [Hz]', 'Amplitude [m/s^2]']
        )
        best_harmonic_series.index += 1
        display(tabulate(best_harmonic_series, headers='keys', tablefmt='html'))
    
        # Plot found harmonic series
        fig, ax = plt.subplots(1, 8, figsize=(30, 4))
        for i in range(8):
            s = h_series[i+1]
            if i == max_harmonic_amp_idx:
                ax[i].xaxis.label.set_color('red')
    
            ax[i].plot(freqs, Pxx.T[i_window])
            ax[i].scatter([x[0] for x in s], [x[1] for x in s], marker='^', color='red')
            ax[i].set_xlabel('Frequency [Hz]')
    
        plt.show()

def show_spectra_largest_amplitudes(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))

        i_window = window_idx(t)
        x_fundamental = get_max_frequency(freqs, Pxx, i_window)
        peaks = get_peaks(freqs, Pxx, i_window)
        
        display(Markdown(f'- *Fundamental frequency:* {x_fundamental} Hz'))
        display(tabulate(peaks.head(5), headers='keys', tablefmt='html'))


def compare_limited_specrograms(spectrograms: list, axis: str, t: float):
    fig, ax = plt.subplots(len(DATASET), 1, figsize=(20, 20), sharey=True)
    i = 0
    for name, ts in DATASET:
        signal = ts[axis].loc[t:t+NFFT/Fs].to_numpy()
        pxx = np.abs(np.fft.rfft(signal) / len(signal)) 
        freqs = np.fft.fftfreq(len(signal), d=1/Fs)[:len(pxx)]
        #ilast = len(freqs[freqs < F_LIMIT])
        
        ax[i].plot(freqs, pxx)
        ax[i].grid(True)
        ax[i].set_xlabel('Frequency [Hz]')
        ax[i].set_ylabel('Amplitude [m/s^2]')
        ax[i].set_xlim(0, F_LIMIT)
        ax[i].set_ylim(0, 0.4)
        ax[i].set_title(name)
        i += 1


def spectrogram_energy_left_cumulative(spectrograms: list, axis: str, t: float):
    fig, ax = plt.subplots(len(DATASET), 1, figsize=(20, 20), sharey=True)
    i = 0
    for name, ts in DATASET:
        signal = ts[axis].loc[t:t+NFFT/Fs].to_numpy()
        pxx = np.abs(np.fft.rfft(signal) / len(signal)) 
        freqs = np.fft.fftfreq(len(signal), d=1/Fs)[:len(pxx)]
        
        ax[i].plot(freqs, np.cumsum(pxx) / np.sum(pxx))
        ax[i].grid(True)
        ax[i].set_xlabel('Frequency [Hz]')
        ax[i].set_ylabel('Cumulative energy [%]')
        #ax[i].set_xlim(0, 10000)
        ax[i].set_title(name)
        i += 1

Compare mafaulda faults

In [ ]:
compare_limited_specrograms(DATASET, 'x', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [ ]:
compare_limited_specrograms(DATASET, 'y', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [ ]:
compare_limited_specrograms(DATASET, 'z', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image

Compare cumulative sums

In [ ]:
spectrogram_energy_left_cumulative(DATASET, 'x', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [ ]:
spectrogram_energy_left_cumulative(DATASET, 'y', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [ ]:
spectrogram_energy_left_cumulative(DATASET, 'z', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image

Spectrogram in X axis

In [ ]:
x_spectra = get_spectrograms(DATASET, 'x')

ch1.tsv (X axis, Fundamental = 48.0315 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

ch2.tsv (X axis, Fundamental = 48.0315 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

ch3.tsv (X axis, Fundamental = 48.0315 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

ch4.tsv (X axis, Fundamental = 48.0315 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

ch5.tsv (X axis, Fundamental = 48.1299 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

Spectrogram detail in X axis

In [ ]:
show_spectrogram_detail(x_spectra, 'x', T_SEC)

ch1.tsv (X axis @ 5s)

No description has been provided for this image

ch2.tsv (X axis @ 5s)

No description has been provided for this image

ch3.tsv (X axis @ 5s)

No description has been provided for this image

ch4.tsv (X axis @ 5s)

No description has been provided for this image

ch5.tsv (X axis @ 5s)

No description has been provided for this image

Peaks in frequency spectrum in X axis

  • MMS peak finder algorithm
In [ ]:
show_mms_peaks(x_spectra, 'x', T_SEC)

ch1.tsv (X axis @ 5s)

No description has been provided for this image

ch2.tsv (X axis @ 5s)

No description has been provided for this image

ch3.tsv (X axis @ 5s)

No description has been provided for this image

ch4.tsv (X axis @ 5s)

No description has been provided for this image

ch5.tsv (X axis @ 5s)

No description has been provided for this image

Harmonic series detection in X axis

In [ ]:
# show_harmonic_series(x_spectra, 'x', T_SEC)
In [ ]:
show_spectra_largest_amplitudes(x_spectra, 'x', T_SEC)

ch1.tsv (X axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.0420568 1
1 61.03520.003247851.25
2114.746 0.001998762.35
3168.457 0.001680383.45
4214.844 0.001594874.4

ch2.tsv (X axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.0376015 1
1 95.21480.003458511.95
2178.223 0.002037243.65
3366.211 0.001917027.5
4 73.24220.001673461.5

ch3.tsv (X axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.8281 0.040525 1
1 2.441410.00322946 0.05
2354.004 0.00226791 7.25
3576.172 0.0019375811.8
4270.996 0.00178902 5.55

ch4.tsv (X axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.0435111 1
1 24.41410.00283807 0.5
2 12.207 0.00192549 0.25
3390.625 0.00181199 8
4590.82 0.0016555912.1

ch5.tsv (X axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.0395062 1
1136.719 0.002658242.8
2380.859 0.001789497.8
3 63.47660.001596881.3
4441.895 0.0014002 9.05

Spectrogram in Y axis

In [ ]:
y_spectra = get_spectrograms(DATASET, 'y')

ch1.tsv (Y axis, Fundamental = 48.0315 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

ch2.tsv (Y axis, Fundamental = 48.0315 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

ch3.tsv (Y axis, Fundamental = 48.0315 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

ch4.tsv (Y axis, Fundamental = 48.0315 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

ch5.tsv (Y axis, Fundamental = 48.1299 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

Spectrogram detail in Y axis

In [ ]:
show_spectrogram_detail(y_spectra, 'y', T_SEC)

ch1.tsv (Y axis @ 5s)

No description has been provided for this image

ch2.tsv (Y axis @ 5s)

No description has been provided for this image

ch3.tsv (Y axis @ 5s)

No description has been provided for this image

ch4.tsv (Y axis @ 5s)

No description has been provided for this image

ch5.tsv (Y axis @ 5s)

No description has been provided for this image

Peaks in frequency spectrum in Y axis

In [ ]:
show_mms_peaks(y_spectra, 'y', T_SEC)

ch1.tsv (Y axis @ 5s)

No description has been provided for this image

ch2.tsv (Y axis @ 5s)

No description has been provided for this image

ch3.tsv (Y axis @ 5s)

No description has been provided for this image

ch4.tsv (Y axis @ 5s)

No description has been provided for this image

ch5.tsv (Y axis @ 5s)

No description has been provided for this image

Harmonic series detection in Y axis

In [ ]:
# show_harmonic_series(y_spectra, 'y', T_SEC)
In [ ]:
show_spectra_largest_amplitudes(y_spectra, 'y', T_SEC)

ch1.tsv (Y axis @ 5s)

  • Fundamental frequency: 144.04296875 Hz
f y 1x
0144.043 0.0321256 1
1 31.73830.003822270.220339
2190.43 0.002845641.32203
3651.855 0.001836564.52542
4122.07 0.001695050.847458

ch2.tsv (Y axis @ 5s)

  • Fundamental frequency: 144.04296875 Hz
f y 1x
0144.043 0.0349968 1
1151.367 0.003043431.05085
2 26.85550.002523310.186441
3 41.50390.002241720.288136
4532.227 0.001896053.69492

ch3.tsv (Y axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.0306644 1
1 58.59380.003184951.2
2480.957 0.002519779.85
3109.863 0.001956012.25
4295.41 0.001684216.05

ch4.tsv (Y axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.0393909 1
1 31.73830.004347590.65
2 68.35940.002855361.4
3 21.97270.002259870.45
4373.535 0.001885197.65

ch5.tsv (Y axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.0379273 1
1432.129 0.00364708 8.85
2576.172 0.0025643511.8
3368.652 0.00188733 7.55
4673.828 0.0015862113.8

Spectrogram in Z axis

In [ ]:
z_spectra = get_spectrograms(DATASET, 'z')

ch1.tsv (Z axis, Fundamental = 0.1969 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

ch2.tsv (Z axis, Fundamental = 0.0984 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

ch3.tsv (Z axis, Fundamental = 6.2992 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

ch4.tsv (Z axis, Fundamental = 0.5906 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

ch5.tsv (Z axis, Fundamental = 0.8858 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

Spectrogram detail in Z axis

In [ ]:
show_spectrogram_detail(z_spectra, 'z', T_SEC)

ch1.tsv (Z axis @ 5s)

No description has been provided for this image

ch2.tsv (Z axis @ 5s)

No description has been provided for this image

ch3.tsv (Z axis @ 5s)

No description has been provided for this image

ch4.tsv (Z axis @ 5s)

No description has been provided for this image

ch5.tsv (Z axis @ 5s)

No description has been provided for this image

Peaks in frequency spectrum in Z axis

In [ ]:
show_mms_peaks(z_spectra, 'z', T_SEC)

ch1.tsv (Z axis @ 5s)

No description has been provided for this image

ch2.tsv (Z axis @ 5s)

No description has been provided for this image

ch3.tsv (Z axis @ 5s)

No description has been provided for this image

ch4.tsv (Z axis @ 5s)

No description has been provided for this image

ch5.tsv (Z axis @ 5s)

No description has been provided for this image

Harmonic series detection in Z axis

In [ ]:
# show_harmonic_series(z_spectra, 'z', T_SEC)
In [ ]:
show_spectra_largest_amplitudes(z_spectra, 'z', T_SEC)

ch1.tsv (Z axis @ 5s)

  • Fundamental frequency: 144.04296875 Hz
f y 1x
0144.043 0.006491461
1 80.56640.0040296 0.559322
2151.367 0.002990131.05085
3205.078 0.001979871.42373
4 70.80080.001741170.491525

ch2.tsv (Z axis @ 5s)

  • Fundamental frequency: 2.44140625 Hz
f y 1x
0 2.441410.00865982 1
1144.043 0.00525229 59
2153.809 0.00282942 63
3 68.3594 0.00195273 28
4268.555 0.00166528 110

ch3.tsv (Z axis @ 5s)

  • Fundamental frequency: 36.62109375 Hz
f y 1x
0 36.62110.009750711
1102.539 0.004403672.8
2 92.77340.002610922.53333
3187.988 0.002023565.13333
4263.672 0.001691047.2

ch4.tsv (Z axis @ 5s)

  • Fundamental frequency: 14.6484375 Hz
f y 1x
0 14.64840.00765251 1
1 21.97270.00464057 1.5
2 85.44920.00267801 5.83333
3183.105 0.0023300212.5
4134.277 0.00171221 9.16667

ch5.tsv (Z axis @ 5s)

  • Fundamental frequency: 17.08984375 Hz
f y 1x
0 17.08980.00767103 1
1 46.38670.00453737 2.71429
2124.512 0.00275187 7.28571
3183.105 0.0017487910.7143
4231.934 0.0014846413.5714

Histogram

In [ ]:
axis = ['x', 'y', 'z']
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts[axis].hist(figsize=(10, 5), grid=True, bins=50)
    plt.show()

ch1.tsv

No description has been provided for this image

ch2.tsv

No description has been provided for this image

ch3.tsv

No description has been provided for this image

ch4.tsv

No description has been provided for this image

ch5.tsv

No description has been provided for this image
In [ ]:
axis = ['x', 'y', 'z']
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts[axis].boxplot(figsize=(10, 5))
    plt.show()

ch1.tsv

No description has been provided for this image

ch2.tsv

No description has been provided for this image

ch3.tsv

No description has been provided for this image

ch4.tsv

No description has been provided for this image

ch5.tsv

No description has been provided for this image

Orbitals of all cross sections

In [ ]:
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    fig, ax = plt.subplots(1, 3, figsize=(20, 4))

    for i, col in enumerate([('x', 'y'), ('x', 'z'), ('y', 'z')]):
        ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
        ax[i].grid(True)
        ax[i].set_xlabel(col[0].upper())
        ax[i].set_ylabel(col[1].upper())
        ax[i].grid(True)
    plt.show()       # plt.savefig('orbitals.png')

ch1.tsv

No description has been provided for this image

ch2.tsv

No description has been provided for this image

ch3.tsv

No description has been provided for this image

ch4.tsv

No description has been provided for this image

ch5.tsv

No description has been provided for this image

Orbitals of 1x harmonic frequency

In [ ]:
x_spectra_by_name = {spec[0]: spec for spec in x_spectra}
y_spectra_by_name = {spec[0]: spec for spec in y_spectra}
z_spectra_by_name = {spec[0]: spec for spec in z_spectra}
t = 5
space = 5

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    fig, ax = plt.subplots(1, 3, figsize=(20, 4))

    name, freqs, Pxx = x_spectra_by_name[name]
    x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
    name, freqs, Pxx = y_spectra_by_name[name]
    y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    name, freqs, Pxx = z_spectra_by_name[name]
    z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    try:
        ts['x_1x'] = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
        ts['y_1x'] = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
        ts['z_1x'] = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
    except ValueError:
        continue
    
    for i, col in enumerate([('x_1x', 'y_1x'), ('x_1x', 'z_1x'), ('y_1x', 'z_1x')]):
        ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
        ax[i].grid(True)
        ax[i].set_xlabel(col[0].upper())
        ax[i].set_ylabel(col[1].upper())
        ax[i].grid(True)
    
    plt.show()       # plt.savefig('orbitals_1x.png')

ch1.tsv

No description has been provided for this image

ch2.tsv

ch3.tsv

No description has been provided for this image
No description has been provided for this image

ch4.tsv

No description has been provided for this image

ch5.tsv

No description has been provided for this image
In [ ]:
t = 5
space = 8

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))

    name, freqs, Pxx = x_spectra_by_name[name]
    x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
    name, freqs, Pxx = y_spectra_by_name[name]
    y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    name, freqs, Pxx = z_spectra_by_name[name]
    z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    try:
        x = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
        y = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
        z = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
    except ValueError:
        continue

    ax = plt.figure().add_subplot(projection='3d')
    ax.scatter(x, y, z, zdir='z', s=1, color='navy')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_xlim(-2, 2)
    ax.set_ylim(-2, 2)
    ax.set_zlim(-2, 2)
    ax.zaxis.labelpad = -0.7
    plt.show()

ch1.tsv

No description has been provided for this image

ch2.tsv

ch3.tsv

No description has been provided for this image

ch4.tsv

No description has been provided for this image

ch5.tsv

No description has been provided for this image